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Abstract 

We present a numerical study on the dynamics of imbibition fronts in porous 
media using a pipe network model. This model quantitatively reproduces the 
anomalous scaling behavior found in imbibition experiments [Phys. Rev. E 
52, 5166 (1995)]. Using simple scaling arguments, we derive a new identity 
among the scaling exponents in agreement with the experimental results. 
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Self-affine interfaces are found in a variety of phenomena such as two-phase flow in 
porous media, thin film deposition, flame fronts, etc. JTJ. In particular, scaling properties 
of imbibition fronts have been the focus in many experimental studies Static scaling 

behavior has been suggested to be described by a directed depinning percolation model 
To explain dynamic properties, one need to go beyond local models and properly deal 
with the effects of long-range coupling through pressure. Much progress has been obtained 
using theoretical approaches including consideration of capillary waves ||, a Flory-type 
scaling argument |§ and very recently a phase-field model fT0|- These studies contributed 



to our better understanding of the scaling behavior of imbibition. However, exponents 
provided in general do not compare satisfactorily with the experimental values. In fact, 
of the four exponents determined experimentally by Horvath and Stanley || with good 
precision, none has previously been explained analytically or numerically, and no existing 
exponent identity applies. In this work, we successfully reproduce all these exponents using 
a pipe network model once the relevant model parameters are properly fine tuned. A new 
exponent identity among these exponents is also presented. In addition, we have found an 
empirical relation between the spatial and temporal correlations. 

Horvath and Stanley studied imbibition fronts climbing up in vertical filter paper sheets 
which move continuously downward into a water container at constant speed v ||. The 
interface height at horizontal coordinate x and time t is denoted by h(x,t). The average 
height h and width w were demonstrated to scale as 

v ~ h and w ~ v~ K (1) 



The temporal height-height correlation function C(t)=< [h(x,t' + t) — h(x,t')} 2 > x / 2 JTT 
was shown to follow a scaling form, 

C(t) = v- K f(tv {et+R)/ P), (2) 

where f(u) ~ u 13 for small u, and it converges to a constant for large u. The exponents were 
found to be 



Q = 1.594 , k = 0.48 , 6 t = 0.37 , and (3 = 0.56. (3) 

The first exponent has particularly strong implications because Darcy's law, which is well 
established for many porous media ||12|| , implies Q = 1. The same result was also found in 
numerous numerical studies | JTTj|rj~3| , [T^ ] . Physically, the anomalous value Q = 1.594 suggests 



that as the interface height increases, its propagation speed becomes slower than simple con- 
siderations of capillary pressure and viscous drag would suggest. This cannot be explained 
by gravity, evaporation, deformation of wetted paper |H|, or inertial effects [12]| since they 



involve intrinsic length or time scales and can only introduce transients or cut-offs to Darcy's 
prediction. 

Our network model consists of cylindrical pipes connecting volumeless nodes on a two- 
dimensional square lattice with periodic boundary conditions in the horizontal direction. 
According to Possiuelle's law |12|], the flow rate in a pipe of radius r is Q = 7rr 4 Ap/8/i<i 
where Ap is the pressure difference between two points separated by d. The fluid viscosity is 
denoted by fi. The pressure at the water level in the tank is maintained at the atmospheric 
value. The capillary pressure is 27 cos <\>jr in each partially filled pipe where 7 and <p are the 
surface tension and the contact angle respectively. Gravity is neglected. This is because in 
Ref. |§, scalings were reported only for h smaller than 65% of the maximum value limited 
by gravity at v = [pL5|| . We have checked numerically that gravity is negligible in this 
range, even though it can lead abruptly to the breakdown of Eq. ([!]) at larger h and finally 
to the pinning transition. We calculate the pressure at every wetted node above the water 
level of the tank by solving Kirchoff's equations using a locally adaptive over-relaxation 
method. The propagation of the interface during the associated adaptive time step is hence 
computed. For computational convenience, possible deviation from Possiuelle's law at the 
junction of the pipes is neglected as usual |L3] and menisci cannot retreat after reaching a 
node. Since in real situation air can escape from either side of the sheet, trapping of air is 
irrelevant and its flow is not simulated. 

The pipe model at this point is similar to those in previous studies [T2|-|H| most of 
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which focused on flow inside porous rocks. Tenuous percolation type wetting patterns are 
obtained. To simulate more compact patterns observed in paper, we therefore implement a 
new wettability rule to impose a strong effective surface tension in our network. Now, water 
can enter an empty pipe from a wetted node only if the neighborhood is sufficiently wet. 
Specifically, of the 8 most nearby nodes (the 4 nearest and the 4 next nearest neighbors) 
forming a loop around the node under consideration, there must exist a connected subset 
of 5 which are already wet. This is the most stringent criterion for the wettability of pipes 
connected to a wetted node without halting the imbibition entirely or involving consideration 
of further neighbors. It leads to a very strong effective surface tension. Different and more 



detailed wettability rules were discussed in Ref. [[L4|]. In addition, fluid pathways in randomly 
arranged fibers have complicated geometry and there are abundance of narrow bottlenecks 
and large pores. We therefore propose another important characteristic of our network, 
namely large spatial fluctuations in the local properties. Its implementation will be discussed 
later. 

The model was first tested in the presence of small fluctuations. Adopt a unit lattice 
spacing and specify the fluid properties by /i = 7cos0 = 1 [Tj|. We assume random pipe 



radii in the range [0,0.5] uniformly distributed JT7|] . This distribution is broader than those 
used by others |l3|Jl4l| . The Darcy's prediction Q = 1 is readily verified. 

To simulate large fluctuations expected in real paper, we simply replace a fraction p w of 
pipes by wider ones with radius r w . For sufficiently large p w and r w , Q > 1 is obtained. Both 
exponents Q and k then depend non-trivially but continuously on p w and r w . Fortunately, we 
have found that nice power-laws with the experimental values of Q and k can be reproduced 



if we put p w = 0.18 and r w = 2.0 [17|. We will adopt these parameters in the remainder of 
this work ||18|| . Snapshots of simulations of imbibition in stationary sheets are shown in Fig. 

We have simulated imbibition in moving paper sheets considered in Ref. || . This involves 
continuous upward shifting of the level on the lattice which represents the contact line of 
the paper sheet with the water in the reservoir. The interface height h and width w are 



measured after steady state has been attained and are averaged over 3 independent runs 
using lattices of width 200 and height 1000. Figures 0(a) and (b) plot respectively v and w 
against h. The linearity observed in log- log plots verifies the scaling relations in Eq. ([[]). 
We get Q — 1.62 ± 0.05 and k = 0.49 ± 0.03 close to the experimental values in Eq. (|3]) due 
to the fine tuning of the pipe distribution mentioned above. The quoted errors represent 
only the uncertainties in the linear fits. 

The spatial correlation C(l)=<[h(x+l, t)-h{x, t)] 2 > 1//2 and the temporal counterpart C(t) 
defined earlier are computed and plotted in Figs. ||](a) and (b) respectively. The initial linear 
regions corresponding to C(l) ~ l a and C(t) ~ t& for small I and t give a = 0.61 ± 0.01 and 
(3 = 0.63±0.01 respectively. This value of j3 is in reasonable agreement with the experimental 
value 0.56. We have verified that C(t) in Fig. [|(b) follows the experimentally motivated 
scaling forms in Eq. @ and we obtain 9 t = 0.38 ± 0.02 in excellent agreement with the 
experimental value 0.37. Similarly, the spatial correlation follows an analogous scaling form 
C(l) = v~ K g{lv^ l+K ^ a ) where g(u) ~ u a for small u and it converges to a constant for large 
u. A further result is obtained by first reparametrizing C(t) by the vertical displacement 
z — vt of the sheet. Equation (0) then becomes C(z) = v" K f(zv^ ez+K ^ 13 ), where 8 Z = 8 t — (3. 
We note that both C(l) and C(z) can be collapsed individually using a common set of 
exponents a = f3 = 0.62 and Q\ = 9 Z = —0.24. Furthermore, they can all be collapsed 
together as shown in Fig. |3](c) verifying the empirical identities 

a — (3 , 0i = 9 Z and f{cu) = g{u) (4) 

where c = 1.6. 

Next, we simulate imbibition in stationary paper sheets. The interface height h and 
width w averaged over 7 lattices of width 1000 are found to scale as 

h~t e , and w~ t Pw (5) 

where e = 0.374 ± 0.005 and (3 W = 0.29 ± 0.01. The interface speed v(t) = dh/dt hence 
computed from direct numerical differentiation on our data is plotted as a function of h in 
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Fig. H](a). The width w is also plotted against h in Fig. 0(b). Direct comparison with the 
moving sheet results is possible because the paper speed v in that case is also the average 
interface speed in the paper frame. These plots show that the scaling relations in Eq. ([!]) 
apply to stationary sheets with the same exponents as well. Therefore, we can derive Eq. 
(H) from Eq. ([5]) and vice versa. In the derivation, one obtains 

and f3 w = — — , (6) 



tt + l ^ + 1 

which describe our exponents accurately. 

For a stationary sheet, w at a given h is only about 80% of the corresponding value 
for a moving sheet as observed from Figs. 0(b). The roughness has therefore not yet 
fully developed for the given height h. In contrast, for the less noisy case the roughness has 
saturated completely and can be determined solely from h [|Hj] . More importantly, our result 
indicates that the roughness gets neither closer to nor further from saturation as h increases 
since w is always 80% of the saturated value. This constancy of the degree of saturation 
implies that the system has a unique dynamic time scale dictating both the growth of the 
roughness and the steady state dynamics of a roughened interface. The time it takes to 
develop a width w is w l ^ w from Eq. (|5|). It is thus proportional to the relaxation time 
extracted from the temporal correlation function u~( e *+ /t )// 3 ~ w { 9 t+n)/(/3K) f rom gq S _ (g) anc [ 
([!]). Therefore, 1/ j3 w = (9 t + k)/{(3k). Applying Eq. (|5|), we have 

/3 = n(0 t + «)/(n + i). (7) 

This identity is particularly important because all values have been measured experimentally 
in Ref. 0. Inserting the experimental values into the r.h.s. of Eq. ([?]), we get (3 = 0.52 
in reasonable agreement with the experimental value 0.56. Using our numerical estimates 
instead, we get (3 = 0.54 which is a little smaller than the numerically found value 0.61. 
In the later case we relate the discrepancy to errors in the numerical determination of (3. 
This may be due to lattice discretization effects which becomes more important at C(t) ^ 1. 
The spatial counterpart a should hence suffer a similar problem. Furthermore, noticeably 
different values of a and f3 are obtained if correlations of higher moments are used [|l(| . 
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The good agreement between our numerical results and the experimental ones strongly 
supports that our model has captured the essential physics in the imbibition process. How- 
ever, some further points are yet to be considered. First, despite the nice power-law fits and 
data collapses of the relevant quantities observed, existence of extraordinarily slow crossover 
effects masking different asymptotic scalings should be cautioned. Second, our model re- 
produces the experimental results after fine tuning of the radius distribution. We have also 
found a few other completely different distributions which also lead to similar exponents 
upon fine tuning fl8| . The exponents nevertheless are not robust with respect to changes in 



the details of the models. Using distributions other than the fine tuned ones, the exponents 



are in general different and sometimes scalings may not even hold fll8| . We expect that 
there is a yet unknown selection mechanism so that models generating the experimental 
exponents are preferred. This is currently under active investigations. However, one cannot 
determine apriori a correct realistic pipe distribution. This is because all pipe networks 
are simplified models of fluid pathways which are indeed very different from those of pa- 
per. Realistic simulations, for example, using lattice-Boltzmann method with sophisticated 
boundary conditions in 3 dimensions unfortunately covers only microscopic regions and 
thus cannot be applied to study the macroscopic scalings. 

Finally, local models of imbibition have been classified into isotropic and anisotropic 
universality classes exemplified respectively by the random field Ising model (RFIM) and 
the directed percolation depinning (DPD) model |J. The anisotropy in the DPD model is 
due to a solid-on-solid condition. In contrast, for our model both the wettability rule and 
the flow dynamics treat the horizontal and vertical directions equivalently. Concerning the 
local symmetry, it is more closely related to the isotropic class. It was found that a = 1 and 
Pw — 3/4 for the isotropic class || while e = a = (3 W ~ 0.633 for the anisotropic one at 
criticality. This is to be compared with e ~ 0.374, a ~ 0.61 and (3 W ~ 0.29 for our model. 
It is easy to see that the much larger values of e and (3 W for DPD and also f3 w for RFIM 
are due to the locality of the interections. We believe that the agreement of a between 
DPD and our model is just a coincidence. The morphologies of the surfaces are indeed 



visually quite different. Previous works on non-local models including pipe networks and 
the phase-field model [|l(],[n|,|rjj] are in a weak fluctuations regime and are all consistent only 
with Darcy's law. Our pipe network with strong fluctuations is the only model exhibiting 
distinctly different scalings in good agreement with the experimental ones in Ref. [Q. We 
therefore suggest that it belongs to a new universality class. 

In conclusion, we have simulated imbibition in paper using a pipe network model and 
reproduced all scaling behaviors observed experimentally in Ref. |§. We obtain 

n = 1.62 , k = 0.49 , t = 0.38 and (3 = 0.63 (8) 

which is to be compared with Eq. ([3D. The model displays rich behaviors, and can be tuned 
to reproduce the experimentally determined values of Q and k. Then (3 and 9 t turn out 
respectively in reasonable and excellent agreement with the experimental values. Assuming 
a single time scale in the dynamics, we have presented a new exponent identity (Eq. (]?])) 
which is justified by the experimental values. Two other exponent identities in Eq. (^) 
relating the moving and stationary paper cases are deduced and verified numerically. Further 
identities in Eq. ([|) for exponents and scaling functions are suggested empirically based on 
a data collapse between the spatial and temporal correlations. 
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FIGURES 




FIG. 1. (a) Snapshot of a small scale simulation showing individual pipes. Wet and dry regions 
are shaded in black and grey respectively. Radii of pipes have been scaled down, (b) Snapshot of 
the wetted region in a larger scale simulation on a lattice of width 500. The average height is 36.8. 
The arrows provided only to demonstrate the ratio of the scales of the two figures. 
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FIG. 2. Log-log plots of (c) v and (b) w against h. In (a), the fitted lines have slopes -1.62 and 
-1.65 for the moving and stationary case respectively. In (b), they are 0.785 and 0.786. 
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FIG. 3. Log-log plots of (a) spatial correlation, (b) temporal correlation and (c) simulta- 
neous data collapse of all correlation curves in both (a) and (b) for paper speed v given by 
v/10- 5 = 2, 3, 4, 6, 9, 14, 20, 30 and 40. 
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